Properties of the reaction front in a reaction-subdiffusion process 
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o: 

f^ ^ We study the reaction front for the process A + B -^ C in which the reagents move subdiffusively. 

fSJ . We propose a fractional reaction-subdiffusion equation in which both the motion and the reaction 

, , ' terms are affected by the subdiffusive character of the process. Scaling solutions to these equations 

O ' are presented and compared with those of a direct numerical integration of the equations. We find 

.^*H , that for reactants whose mean square displacement varies sublinearly with time as (r^) ~ f^, the 

^ ' scaling behaviors of the reaction front can be recovered from those of the corresponding diffusive 

0^ ' problem with the substitution t -^ t'' . 

'T? ; I. INTRODUCTION 

o ■ 

^5 I It is very well known that diffusion-limited binary reactions in low dimensions exhibit "anomalous" rate laws, that 

^ ■ is, rate laws that deviate from the quadratic forms that usually characterize binary reactions. It is also known that 
I '. ' these anomalies occur because diffusion is not an effective mixing mechanism, especially in constrained geometries, 
C^ . whereas the quadratic rate laws assume perfect mixing. The shortcomings of diffusion as a mixing process lead to 
c/3 ' the spontaneous occurrence of spatial order and spatial structures in terms of which the anomalous rate laws are well 
^_i understood [l|. Broadly speaking, spatial order leads to a slowing down of the reaction because it tends to decrease 
2 or deplete the interfacial areas where reactions can occur. 

r^ The anomalies exhibit themselves in different ways according to the design of the experiment, and the critical 

' , dimension above which the behavior of the reaction is "normal" also varies from one experiment to another. Thus, 

pi ^ for example, in an irreversible batch reaction A + B ^ C, the decay of reactant concentrations deviates from (is 
Q i slower than) the mean field behavior p ^ t^^ up to dimension d = 4. This anomaly is accompanied by and directly 
^ ' associated with a sharp segregation of species as the majority species in any region eliminates the minority species and 
'"^ , diffusion is too slow to replenish the minority population effectively. If there are random sources of reactants so that 

. • the system achieves a steady state, species segregation is observed in the steady state up to dimension d — 2 [llQl- 
^ One of the more accessible experimental setups to observe these anomalies is that of a reaction front in which 

^\ initially one reactant uniformly occupies all the space to one side of a sharp front, and the other reactant occupies the 
C^ other |3]- As the reaction proceeds, one can monitor the concentration profiles of both reactants and of the product. 
■^ Their evolution reflect the anomalies. For the reaction-diffusion case the front evolution has been extensively studied 

Jtj theoretically and confirmed experimentally |j|. 

Our interest lies in understanding binary reaction kinetics in media where the reactants move subdiffusively. Sub- 
diffusive motion is particularly important in the context of complex systems such as glassy and disordered materials, 
in which pathways are constrained for geometric or for energetic reasons. It is also particularly germane to the way in 
j^ ' which experiments in low dimensions have to be carried out. Such experiments must avoid any active or convective 
Cj or advective mixing so as to ensure that any mixing is only a consequence of diffusion. To accomplish this usually 

I . requires the use of gel substrates and/or highly constrained geometries (the first gel-free experiments were carried out 
'^ ' recently |a, Q)- Under these circumstances it is not clear whether the motion of the species is actually diffusive, or if 
C , it is in fact subdiffusive. Indeed, a recent detailed discussion on ways to extract accurate parameters and exponents 
from such experiments concludes that at least the experiments presented in that work, carried out in a gel, refiect 
subdiffusive rather than diffusive motion |7(. 

Subdiffusive motion is characterized by a mean square displacement that varies sublinearly with time. 
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with < 7 < 1. For ordinary diffusion 7 = 1, and Ki = D is the ordinary diffusion coefficient. Subdiffusion 
is not modeled in a universal way in the literature. Among the more successful approaches to the subdiffusion 
problem have been continuous time random walks with non-Poissonian waiting time distributions ja, 0, llfl HHi 
and fractional dynamics approaches in which the ordinary diffusion operator is replaced by a generalized fractional 
diffusion operator 10, \\2i Il3 . \\M. \W{ . The relation between the two has also been discussed. In particular, the 
fractional dynamics formulation that leads to the mean square displacement Q can be associated with a continuous 
time random walk with a waiting time distribution between steps which at long times behaves as 

4,{t)^t-^-\ (2) 



In our work we adopt the fractional dynamics approach. 

Since subdiffusion is an even poorer mixing mechanism than ordinary diffusion, we expect that the diffusive anoma- 
hes observed in constrained geometries will in some sense be exacerbated in a subdiffusive environment. However, the 
outcome is not entirely clear or predictable. On the one hand, slower motion naturally leads to a slowing down of the 
reaction. On the other, slower motion also delays the elimination of the minority species in a region, and therefore 
may reduce the spatial segregation effects observed in the diffusive case. 

In this paper we focus on the time dependence of the front geometry, and, in particular, on the exponents that 
describe the evolution of a number of quantities that define the reaction front. We expect there to be anomalies that 
differ from those of the diffusive problem. In Sec. |n] we formulate our model and present our scaling solutions. In 
Sec, mil we compare our scaling results with those obtained by direct numerical integration of our model equations. 
We end with some comments on our results in Sec. IIVI 



II. EVOLUTION OF SUBDIFFUSIVE FRONT 

We start with a system of A particles on one side and B particles on the other of a sharp linear front, defined to 
lie perpendicular to the x axis. The particles diffuse and react with a given probability upon encounter. A standard 
mean-field model for the evolution of the concentrations a{x, t) and h{x, t) of A and B particles along x is given by 
the reaction-diffusion equations 

d d^ 

-—a(x,t) = D——:-a(x,t) — kaix,t)h{x,t) 

(3) 



dt ' dx 

d d 

-^b{x,t) = D-—b{x,t)-ka{x,t)b{x,t), 



where D is the diffusion coefficient assumed to be equal for the two species. The initial conditions are that a{x, t) = 
const = ao for x < and a(x, t) — for a; > 0. Similarly, b{x, t) = const — bo for x > and b(x, t) = for a; < 0. With 
these conditions, no matter the dimensionality of the system, the system of equations is effectively one-dimensional. 

The front problem was first analyzed via a scaling description [3J and later refined by a large number of authors 
using more rigorous theoretical and careful numerical approaches (see references in Yuste et al. |j]). One upshot of the 
extensive work is that d = 2 is a critical dimension for the mean field description to be appropriate. Below d = 2 one 
must take into account fiuctuations, neglected in this description, that completely change the outcome of the analysis. 
While it has been assumed that the mean field model holds for the critical dimensions d = 2, Krapivsky ,16] finds 
logarithmic corrections that have also been observed in simulations jl7| . Our program is to generalize the mean field 
description to the subdiffusive case and to apply a scaling description to the problem. One can show that d = 2 is also 
the critical dimension for the subdiffusive problem. Our work is therefore not appropriate for systems of dimension 
lower than two (we do not consider logarithmic corrections). 

In order to generalize the reaction-diffusion problem to reaction-subdiffusion, we must deal with the subdiffusive 
motion of the particles (generalization of the first term in Eq. (Q) and with their reaction rate law (second term). 
We discuss each separately. 

In the fractional diffusion approach to subdiffusive motion one replaces the diffusion operator by the Riemann- 
Liouville operator: 

D-^a{x, t) ^ K, oDJ-'-^aix, t), (4) 

where Kj is the generalized diffusion coefficient that appears in Eq. JQl, and qD^ ^'^ is the Riemann-Liouville operator, 

^D]-V{t)^-^Afdr-I^. (5) 



The reaction-diffusion equations Q are thus replaced by 

— a{x,t) = K^ oDl^'^—a{x,t) - R^{x,t) 



,2 



— 6(x,t) = K^ oDl-''-^bix,t) - R^{x,t). 



(6) 



The reaction term R^{x,t) will be discussed subsequently because certain aspects of the problem are independent of 
the specific form of this term. 



A. Scaling independent of reaction term 

As the reaction proceeds, a depletion zone develops around the front. This is the region where the concentrations 
of reactants are significantly smaller than their initial values. How the width W evolves with time is one of the 
measures typically used to characterize the process. Within this depletion zone lies the so-called reaction zone, the 
region where the concentration c{x, t) of the product C is appreciable. This concentration profile has a width w 
whose variation with time is another characteristic of the evolving reaction. The evolution of the production rate of 
C (which determines the height of the profile of c{x,t) in the reaction zone) is a third measure of the process. To 
find these time dependences we adapt the original scaling approach [j, [l3 to the subdiffusive case, and assume the 
scaling forms 

a{x,t) = t-^aax-Xf)t-") 
b{x,t) = t-'>b{{x~Xf)t-'^) 

for the concentrations and 

R^{x,t)^t-''R^{{x-Xf)t-°') (8) 

for the reaction term. The exponents 6, a, and /i are to be determined from three relations. The scaling forms are 
only valid for x <^ W, that is, well within the depletion zone. 

Two of the three relations needed to fix the scaling exponents do not require further specification of the reaction 
term. Since the reaction zone increases more slowly than the width of the depletion zone (an assumption that ex post 
turns out to be correct), we can focus on the concentration difference u{x,t) — a{x,t) — b{x,t) to deduce the width 
of the latter. The reaction term drops out when one subtracts the equations in Eq. (|SJ), and its form therefore does 
not matter at this point. Generalizing the procedure of Galfi Q, one can scale the resulting equation by measuring 
concentrations in units of ao, time in units of r = l/(A;ao), and length in units of I = (K^t^Y''^, so that the equation 
is simply 

-u{x,t)= oDl-'g^uix,t) (9) 

and the only control parameter is q — bo/aQ in the initial condition: 

uix,0) == 1 for a; < , , 

(10) 
u{x,0) = -g for a; > 0. 



The solution is 



U(x, t) = -q-\ — _Hj J 



t^/2 



(0,1) 



(11) 



2 
where H^ \ is the Fox H- function [l^ [l3| . When 7 = 1 this reduces to the diffusion result Q 

^(^,t)--'Z+^erfc(^). (12) 

From Eq. Hll|l we see that the width of the depletion zone scales as 

W ~ i'^/^ (13) 

i.e., da{x,t)/dx ~ db{x,t)/dx ~ t^''^^. Then, from Eq. IZj), the following relation between scaling exponents follows 
immediately: 

e + a^l (14) 

The second relation follows from the fact that the concentration gradient of A and B leads to a flux of particles 
toward the reaction region. The assumption that the reaction is fed by these particle currents then leads to the 
quasistationary form in the reaction zone 

= K^ oDl-^^^a{x,t) - R^{x,t) 

% (15) 

= if ^ ^D\-''^{x,t) - R^{x,t) , 



which requires that 

^i = + 2a + 1 - 7 . (16) 

For the width of the reaction zone to grow more slowly than the depletion zone caused by subdiffusion requires 
that 

a < 7/2 . (17) 

On the other hand, the quasistationarity condition requires that 

which again leads to Eq. p7|l . 

An experimentally accessible quantity that is independent of R-^ is the location Xf of the point at which the 
production rate of C is largest. This should occur where a{x,t) ^ b{x,t), that is, u{xf,t) ^ 0. The time dependence 
of this equimolar point is found from Eq. Hll() to be 

Xf{t)^Kft'''^ (19) 

where Kf is determined from the equation 
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(0,1) 



(20) 



B. Choice of reaction term and resultant scaling 

Further relations involving the scaling exponents aimed at their expression in terms of model quantities require 
specification of the reaction term. There is a varied literature on this subject, based on a number of different 
assumptions [l9ll2(l l2lL 122. 123 . l24| . Most do not associate a memory with the reaction term. Some assume that, as in 
the case of ordinary diffusion, reactions can simply be modeled by a space-dependent form of the law of mass action, 
e.g., by setting R — ka{x,t)b{x,t). Some of these assumptions may be appropriate if the reaction is very rapid, but 
not if many encounters between reactants are required for the reaction to occur. 

We adopt the viewpoint put forth in a recent theory developed for geminate recombination [23, |2J| but, as the 
authors themselves point out, much more broadly applicable. This theory goes back to the continuous time random 
walk picture from which the fractional diffusion equation can be obtained, and considers both the motion and the 
reaction in this framework. In the context of geminate recombination the authors define a reaction zone and argue 
that a geminate pair within the reaction zone will not necessarily react for any finite intrinsic reaction rate (which 
they call jrc) because one of the particles may leave the zone before a reaction takes place. The dynamics of leaving 
the reaction zone is ruled by the waiting time distribution ipoutit) = ■ip{t)e~^''''^ where ip{t) is the waiting time that 
regulates the rest of the dynamics [cf . Eq. jSJl] , and therefore the reaction rate will acquire a memory that arises 
from the same source as the memory associated with the subdiffusive motion. In the continuum limit this model then 
leads to a reaction-subdiffusion equation in which both contributions have a memory. Seki et al. 23, 24] obtain a 
subdiffusion-reaction equation which at long times corresponds to choosing a reaction term of the form 

R^{x, t) = k aDl"'a{x, t)b{x, t). (21) 

Here "long times" set in very quickly if the reaction zone is narrow and the intrinsic reaction rate small. As noted 
earlier, although the derivation is specifically for geminate recombination, the arguments can be generalized. 
Our full reaction-subdiffusion starting equations on which the remainder of this paper is based then are 



Tra\x,t) — lia\x,t)b[x,t} > 
'^ I 

(22) 



—-a{x,t) = qDI '^ ^K^-^a{x,t) - ka{x,t)b{x,t) 



32 



—b{x,t) ^ oDI^'' ■lK^-—^b{x,t) ka{x,t)b{x,t) 

From the specific reaction term given in Eq. (|22|l we can now obtain the third relation between the scaling exponents 
by balancing the terms within the brackets: 

fj.^2e + l-j. (23) 




FIG. 1: Log-log plots of the width of the product profile iu(i) and the width of the depletion zone W{t) vs t. The squares are 
the numerical results for w{t) and the comparison line has slope 7/6. The circles are the numerical results for W{t) and the 
line has slope 7/2. 



Simultaneous solution of Eqs. 114|) . H16|) . and H23|) finally yields 
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(24) 



III. COMPARISONS WITH NUMERICAL SOLUTIONS 



Elsewhere we have carried out exhaustive numerical simulations of the reaction-subdiffusion system and compared 
our scaling results to those of the simulations for a variety of front properties |3| ■ The simulations were described in 
detail and are particularly subtle when motions are slow and concentrations are low. The agreement with our scaling 
results was in some cases quantitative and in others, particularly where very small numbers are involved, at least 
qualitative. 

Here we present a different comparison, namely, between our exact results (where we have them) and our scaling 
results, and those obtained by direct numerical integration of the reaction-subdiffusion model 25] . This analysis 
provides an insight not only into the accuracy of the assumptions that enter the scaling analysis, but also into the 
length of time (rather short, as it turns out) before the asymptotic scaling results are valid. All of our numerical 
results here have been obtained for the exponent 7 = 0.5 and with the parameters K1/2 — k ^ I. 

Figure n shows the width of the product profile, w{t), predicted to grow as V^^, and the width of the depletion 
zone, W{t), which should grow as f^^ according to our scaling theory. The asymptotic prediction for the former is 
very good after about t — 100 and for the latter after about t — 200. Note that whereas the scaling prediction for 
W{t) is independent of the form of the reaction term, that of the product profile width K;(i) is not. 

The product a(0, t) x 6(0, t) is shown in Fig. |2| The agreement between the numerical integration and the scaling 
result is excellent from the earliest times indicated. 

Figure|31contains a more elaborate set of results showing the space dependence of various quantities at two different 
times. The triangles pointing upward represent the reactant concentration a(x, t) as a function of position x. As 
expected, this profile is large on the left side of the graph (since reactant A was originally located to the left of the 
front). The upper curve (solid triangles) corresponds to time t — 1000 and the lower curve (white triangles), when 
more of A has already reacted, corresponds to t — 5000. On the other side and in complete mirror image (because 
we took both initial concentrations to be equal) are the corresponding results for reactant B. The triangles pointing 
downward represent b{x, t) at the two times, with the same time progression from more reactant to less reactant as time 
moves on. The quantity that we can directly compare to our theory is the difference profile u{x, t) = a(x, t) — h{x, t). 
The numerical integration results are represented by the circles. Way to the left of the diagram u{x,t) p^ a{x,t) and 
way to the right u{x,t) « b{x,t), but near the origin both concentrations contribute to u{x,t). The outer solid curve 




FIG. 2: Log-log plots of the product a(0, t) x b(0, t) vs t. The circles are the numerical results and the line has the slope 27/3 
predicted for this product by our scaling theory. 
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FIG. 3; Profiles a(x,t), b{x,t), and the difference profile u{x,t) = a{x,t) — b{x,t) vs x for two different times. The detailed 
description of this figure is in the text. The theoretical vs numerical comparison is between the circles and the solid curves. 



(through solid circles) is the exact result Eq. Hll|l at t = 1000 and the inner solid curve is for t — 5000. The agreement 
is perfect, as it should be. 

Finally, in Fig. ^ we present results for the time dependence of two quantities frequently used to characterize the 
point of maximum rate of product formation. The point one often sees chosen for this characterization is x — Xf 
such that a{xf,t) = b{xf,t). This time-dependent equimolar point can be found exactly from the solution Eq. I|ll|l . 
Xf = Kft^/"^, where Kf is determined from the equation 
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The circles in the figure are the numerical results for Xf{t) with q — 1/2 and the line is our theoretical prediction with 




FIG. 4: Equimolar point Xf{t) vs t: circles are numerical integration results and the solid line our exact theoretical values. 
The squares are the numerical results for the point of maximum product formation, Xf(t). The difference between the two is 
represented by the diamonds, and the dashed line is our scaling prediction of slope 7/6. 

Kf = 0.482742 .... A second point (and one that is actually a more accurate characterization of the maximum rate 
of product formation) is Xf such that the product a{x f , t)b{x f , t) which determines the reaction rate is a maximum. 
This quantity is more difficult to predict, but we do know that both Xf and Xf are within the reaction zone, which 
scales as xjVl^ . It therefore follows that the distance between these two points should grow as xj — Xf ^ f/^ . But 
since xj ~ P'^, we conclude that Xf ~ Xf only for times for which their difference is negligible, i.e. only for times 
such that i~^^/^ :^ 1. In Fig. 0| the squares are the numerical results for Xf{t), the diamonds represent the difference 
Xf — Xf, and the dashed line has slope 7/6. For the two measures to approach closely one has to go to considerably 
longer times. 

IV. END NOTES 

We have presented a scaling theory for the evolution of a reaction-subdiffusion front for the annihilation reaction 
A + B ^ C. The theory is based on a set of fractional equations in which the diffusion operator is replaced by a 
Riemann-Liouville operator, and the reaction term is adjusted as well to take into account the slower motion of the 
reactants. Elsewhere we presented a comparison of our theoretical results with numerical simulations of the process 4]. 
Here we have presented a detailed comparison of our results with those of numerical integrations of the fractional 
equations [25l | . The familiar rate anomalies of the reaction-diffusion problem are now more pronounced by the slower 
motion of the reactants, with the result that scaling behaviors are recovered from those of the corresponding diffusive 
problem with the substitution t —>■ f^ . 
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